function log_like = LogLikelihood(param)

%% GLOBAL

global Y T1 T2 SM SR Draws M Xc Gmax Gmin or orm ll

%% PARAM def 
cm = param(1);
cr = param(2);
phir = param(3);
c1 = param(4);
phi1 = param(5);
sigem = param(6);
siger = param(7);
sigtheta = param(8);
c2 = param(9);
phi2 = param(10);
cy = param(11);
phim = param(12);
gamma1 = param(13);
gamma2 = param(14);
by = param(15);
b1 = param(16);
b2 = param(17);
cgpa = param(18);
phigpa = param(19);
sigegpa = param(20);

c1or = param(21);
phi1or = param(22);
c2or = param(23);
phi2or = param(24);
b1or = param(25);
b2or = param(26);

% Probit restrictions
sige1 = 1;
sige2 = 1;
sigy = 1;

%% Define individual likelihood

theta = Draws*sigtheta;

Lsm = normalpdf(repmat(SM,1,M)-cm-phim*theta,0,sigem);
Lsr = normalpdf(repmat(SR,1,M)-cr-phir*theta,0,siger);
Lt1 = (normcdf(c1+c1or*or+phi1*theta+phi1or*or.*theta+b1*repmat(Xc,1,M)+b1or*or.*repmat(Xc,1,M), 0, sige1).^repmat(T1,1,M)).*((1-normcdf(c1+c1or*or+phi1*theta+phi1or*or.*theta+b1*repmat(Xc,1,M)+b1or*or.*repmat(Xc,1,M), 0,sige1)).^(1-repmat(T1,1,M)));
Lt2 = (normcdf(c2+c2or*orm+phi2*theta+phi2or*orm.*theta+b2*repmat(Xc,1,M)+b2or*orm.*repmat(Xc,1,M), 0, sige2).^repmat(T2,1,M)).*((1-normcdf(c2+c2or*orm+phi2*theta+phi2or*orm.*theta+b2*repmat(Xc,1,M)+b2or*orm.*repmat(Xc,1,M), 0,sige2)).^(1-repmat(T2,1,M)));
Lgpa = normcdf(cgpa+phigpa*theta-Gmax,0,sigegpa).*(repmat(Xc,1,M)>=Gmax)+(repmat(Xc,1,M)<Gmax & repmat(Xc,1,M)>Gmin).*normalpdf(repmat(Xc,1,M)-cgpa-phigpa*theta,0,sigegpa)+(repmat(Xc,1,M)==Gmin).*(1-normcdf(cgpa+phigpa*theta-Gmin,0,sigegpa));

%% Define Bias

yind = normcdf(cy+theta+by*repmat(Xc,1,M));
yind = yind+0;
t1pr = normcdf(c1+c1or*or+phi1*theta+phi1or*or.*theta+b1*repmat(Xc,1,M)+b1or*or.*repmat(Xc,1,M));
t2pr = normcdf(c2+c2or*orm+phi2*theta+phi2or*orm.*theta+b2*repmat(Xc,1,M)+b2or*orm.*repmat(Xc,1,M));
bias1 = repmat(T1,1,M)-yind;
bias2 = repmat(T2,1,M)-yind;

Ly = (normcdf(cy+theta+gamma1*bias1+gamma2*bias2+by*repmat(Xc,1,M), 0, sigy).^repmat(Y,1,M)).*((1-normcdf(cy+theta+gamma1*bias1+gamma2*bias2+by*repmat(Xc,1,M), 0, sigy)).^(1-repmat(Y,1,M)));

%% Individual Likelihood Contribution
Lm=((((Lsm.*Lsr).*Lt1).*Lt2).*Ly).*Lgpa;

%% log likelihood
ll = (log(mean(Lm,2)));
log_like = -sum(log(mean(Lm,2)));
